######################################################## 
##  Program Name  :/projects/bsi/bioinf_int/s106381.borawork/beagle_test/rpgm/compare_bglprob_geno.R
##  Study Title   : BORA
##  Programmer    : MEM
##  Lead :          de Andrade/Sicotte
##  Creation eDate: Monday, 21 February 2011 04:44 PM CST
##  Function      : Compare BEAGLE imputed genotype probabilities to actual genotype
######################################################## 
indir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/1000samples"

##Get gprobs
setwd(indir)

##Get whole chrom run
whole <- read.table("whole_default/beagle_out_chr22.geno.bgl.gz.gprobs", header=T, as.is=T)
whole[,seq(4,ncol(whole),by=3)] <- 2*whole[,seq(4,ncol(whole),by=3)]+whole[,seq(5,ncol(whole),by=3)]
whole <- whole[,c(1:3, seq(4,ncol(whole),by=3))]

# preparing map file
map <- read.table("./whole_default/markers.txt", sep=" ")
names(map)<-c("rsid","bp","alleleA","alleleB")
map$chr <- 22
rownames(map)<-map[,1]
r2 <- read.table("whole_default/beagle_out_chr22.geno.bgl.gz.r2", header=F, as.is=T)
rownames(r2)<-r2[,1]
map$r2<-r2[rownames(map),2]
map<-map[,c("rsid","bp","chr","r2","alleleA","alleleB")]

#types of edges & windows
gprob <-c("./1mb/whole_gprobs","./2mb/whole_gprobs","./4mb/whole_gprobs","./halfmb/whole_gprobs","./halfquatermb/whole_gprobs","./quatermb/whole_gprobs","./quaterquatermb/whole_gprobs")

for (g in gprob) {
  print(g)
  sp <- gsub(".*[[:punct:]]", "", gsub("mb.*", "mb", g))
  temp <- read.table(g, header=T, as.is=T)
  temp[,seq(4,ncol(temp),by=3)] <- 2*temp[,seq(4,ncol(temp),by=3)]+temp[,seq(5,ncol(temp),by=3)]
  temp <- temp[,c(1:3, seq(4,ncol(temp),by=3))]
  
  ##Get difference
  diff <- (whole[,-(1:3)]-temp[,-(1:3)])

  map[,paste("sumdiff",sp,sep=".")] <- apply(abs(diff),1,sum)
  map[,paste("meandiff",sp,sep=".")] <- apply(abs(diff),1,mean)
  map[,paste("meddiff",sp,sep=".")] <- apply(abs(diff),1,median)
  map[,paste("per1diff",sp,sep=".")] <- apply(abs(diff)>.1,1,mean)
  map[,paste("per2diff",sp,sep=".")] <- apply(abs(diff)>.2,1,mean)
  map[,paste("per5diff",sp,sep=".")] <- apply(abs(diff)>.5,1,mean)
  
}

#for second run with random seed
g<-"whole_random2/beagle_out_chr22.geno.bgl.gz.gprobs"
sp <- "seed"
temp <- read.table(g, header=T, as.is=T)
temp[,seq(4,ncol(temp),by=3)] <- 2*temp[,seq(4,ncol(temp),by=3)]+temp[,seq(5,ncol(temp),by=3)]
temp <- temp[,c(1:3, seq(4,ncol(temp),by=3))]
##Get difference
diff <- (whole[,-(1:3)]-temp[,-(1:3)])
map[,paste("sumdiff",sp,sep=".")] <- apply(abs(diff),1,sum)
map[,paste("meandiff",sp,sep=".")] <- apply(abs(diff),1,mean)
map[,paste("meddiff",sp,sep=".")] <- apply(abs(diff),1,median)
map[,paste("per1diff",sp,sep=".")] <- apply(abs(diff)>.1,1,mean)
map[,paste("per2diff",sp,sep=".")] <- apply(abs(diff)>.2,1,mean)
map[,paste("per5diff",sp,sep=".")] <- apply(abs(diff)>.5,1,mean)
temp<-c()
whole<-c()

##Get breaks
segs <- gsub(".*[[:punct:]]", "", gsub("mb.*", "mb", gprob))
setwd(indir)
segdat <- data.frame("break"=1:3,stringsAsFactors=F)

#for chr 22 only you need to change value for other chrs
tot_num_window<-seq(1,4)

for (seg in segs) 
{
	window_summ<-c()	   	
	for(num in tot_num_window)
	{
		file<-paste(seg,paste("/1/",num,"/markers_segment",num,'.txt',sep=""),sep="")
		print(file)
		m<-read.table(file, header=F, as.is=T)
		window_summ<-rbind(window_summ,c(num,summary(m[,2])[1],summary(m[,2])[6]))
	}
	names(window_summ)<-c("num","min","max")
	midseg<-(window_summ[-1,2]+window_summ[-nrow(window_summ),3])/2
  	segdat[,paste("seg", seg, sep=".")] <- midseg
}



outdir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/result_new_paper_graph_11_01_2011"
setwd(outdir)
#map <- map[!is.na(map$bp),]
png(file="chr22_median_dosage_diff_1000samples_new.png",res=100, width=1000, height=800)
#yrange <- range(c(map$meddiff.1mb, map$meddiff.2mb, map$meddiff.halfmb, map$meddiff.seed),na.rm=T)
yrange<-c(0,0.20)
plot(lowess(map$bp, map$meddiff.seed,f=0.00001), type='l', ylim=yrange, main="Median of Absolute Value of Dosage Difference", xlab="BP", ylab="Median Absolute Difference")
lines(lowess(map$bp, map$meddiff.halfmb,f=0.00001),lty=4, col="green")
lines(lowess(map$bp, map$meddiff.1mb,f=0.00001),lty=2, col="red")
lines(lowess(map$bp, map$meddiff.2mb,f=0.00001), lty=3, col="blue")
lines(lowess(map$bp, map$meddiff.quatermb,f=0.00001), lty=5, col="cyan")
lines(lowess(map$bp, map$meddiff.halfquatermb,f=0.00001), lty=6, col="purple")
lines(lowess(map$bp, map$meddiff.quaterquatermb,f=0.00001), lty=7, col="magenta")

lines(lowess(map$bp[!is.na(map$r2)], yrange[2]-yrange[2]*map$r2[!is.na(map$r2)],f=0.00001), lty=5, col="orange")
abline(v=segdat$seg.1mb, col="red", lty=2)
abline(v=segdat$seg.2mb, col="blue", lty=3)
abline(v=segdat$seg.halfmb, col="green", lty=4)
abline(v=segdat$seg.quatermb, col="cyan", lty=5)
abline(v=segdat$seg.halfquatermb, col="purple", lty=6)
abline(v=segdat$seg.quaterquatermb, col="magenta", lty=7)

#legend("topleft", col=c("black", "green", "red", "blue","cyan","purple","magenta"), lty=c(1,4,2,3,5,6,7), c("Different Seed", "Half MB", "1 MB", "2 MB","QUAT MB","HF_QUAT MB","QU_QUAT MB"), bty='n')
legend("topleft", col=c("orange","black", "green", "red", "blue","cyan","purple","magenta"), lty=c(1,4,2,3,5,6,7), c("Scaled 1-r2, Peaks are\nareas of low quality","Different Seed", "Half MB", "1 MB", "2 MB","QUAT MB","HF_QUAT MB","QU_QUAT MB"), bty='n')
dev.off()
